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Abstract 

Constraints from analyticity are combined with experimental electron-proton scattering 
data to determine the proton charge radius. In contrast to previous determinations, we 
provide a systematic procedure for analyzing arbitrary data without model-dependent 
assumptions on the form factor shape. We also investigate the impact of including 
electron-neutron scattering data, and irir — > NN data. Using representative datasets we 
find r p E = 0.870 ± 0.023 ± 0.012 fm using just proton scattering data; r p E = 0.880±g;g^ ± 
0.007 fm adding neutron data; and r p E = 0.871 ± 0.009 ± 0.002 ± 0.002 fm adding tttt 
data. The analysis can be readily extended to other nucleon form factors and derived 
observables. 



1 Introduction 



The electromagnetic form factors of the nucleon provide basic inputs to precision tests of 
the Standard Model. In particular, the root mean square (RMS) proton charge radius as 
determined by the form factor slopqj , 



1 + + 



(1) 



is an essential input to hydrogenic bound state calculations [HE]- Recent experimental results 
suggest a discrepancy between the charge radius inferred from the Lamb shift in muonic 
hydrogen [3], r% = {r 2 ) p E = 0.84184(67) fm, and the CODATA value, r\ = 0.8768(69) fm, 
extracted mainly from (electronic) hydrogen spectroscopy [I]. The charge radius can also be 
extracted from elastic electron-proton scattering data. The 2010 edition of the Review of 
Particle Physics lists 12 such determinations that span the range of 0.8-0.9 fm [5J, most with 
quoted uncertainties of 0.01-0.02 fm. These determinations correspond to analyses of different 
datasets and different functional forms of G p E (q 2 ) that were fit to the data over a period of 50 
years. 

Extraction of the proton charge radius from scattering data is complicated by the unknown 
functional behavior of the form factor. We are faced with the tradeoff between introducing 
too many parameters (which limits predictive power) and too few parameters (which biases 
the fits). Here we describe a procedure that provides mo del- independent constraints on the 
functional behavior of the form factor. The constraints make use of the known analytic 
properties of the form factor, viewed as a function of the complex variable t = q 2 = —Q 2 . 
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Figure 1: Conformal mapping of the cut plane to the unit circle. 



As illustrated in figure [TJ the form factor is analytic outside of a cut at timelike values 
of t, [6] beginning at the two-pion production threshold, t > 4m^.J§ In a restricted region 
of physical kinematics accessed experimentally, — Q max < t < 0, the distance to singularities 
implies the existence of a small expansion parameter. We begin by performing a conformal 

1 G P E is defined in Section I3TT1 

2 Here and throughout, — 140 MeV denotes the charged pion mass, and mjy = 940 MeV is the nucleon 
mass. 
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mapping of the domain of analyticity onto the unit circle: 

V^cut 



Z\t, t C ut: to) 



t-y/t 



cut 



*0 



y/tcut — t + y/t cut — to 



(2) 



where for this case t. 



cut 



onto z — 0. By the choice tg pt = *cut _ \A 



4m 2 , and to is a free parameter representing the point mapping 



minimized: \z\ < |z| 



2 

max 



= [(i + QLJt c 

0.05 GeV 2 , 0.5 GeV 2 , we find \z\ r _ 



Qmax/^cut J , the maximum value of \z\ is 

- 1]/[(1 + QLxAcut)^ + !]■ For example, with 
0.062, 0.25. Expanding the form factor as 



oo 

E 

fc=0 



2\k 



a k z(q ) 



(3) 



we find that the impact of higher order terms are suppressed by powers of this small param- 
eteJE As we will see below, the coefficients multiplying z k are bounded in size, guaranteeing 
that a finite number of parameters are necessary to describe the form factor with a given pre- 
cision. Figure [2] illustrates the manifestation of this fact in the form factor data. As expected, 
the curvature is smaller in the z variable than in the Q 2 variable. 
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Figure 2: Form factor as a function of Q 2 and as a function of z. Here we choose to = in 
the definition of z, and plot data from [7] for < Q 2 < 0.5 GeV 2 . 



Expansions of the form (|2J) are a standard tool in analyzing meson transition form factors [SI 
El EDI IHl EH HH EES EH CEO HSJ . A complicating feature in the present application to nucleon 
form factors is the contribution of the subthreshold region 4m 2 < t < 4m 2 N in the relevant 
dispersion integral. 

The rest of the paper is structured as follows. In Section [2] we demonstrate the application 
of the z expansion in some illustrative fits and compare it to other expansions that appear in 



3 Physical observables are independent of the choice of to , which can be viewed as the choice of an expansion 
"scheme". |z| max defined in this way gives a convenient estimation of the impact of higher-order terms. 
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the literature. One of the main advantages of the z expansion is that the expansion coefficients 
can be bounded using knowledge about ImG^ in the time-like region. In Section |3] we discuss 
these bounds. In Section H] we discuss several possibilities of reducing the error on the charge 
radius by including more experimental data, namely: high Q 2 data, neutron scattering data, 
and im data. Finally, we discuss our results in Section |5j 



2 Illustrative fits 

Let us consider the six datasets tabulated by Rosenfelder [19] (denoted in [19] as SI, S2, R, 
Bl, B2, M) This will allow us to compare in detail the results of our fit to previous analyses. 
For definiteness, we take all data points in [19] with corrections from magnetic form factor 
contributions A mag < 0.15. The resulting dataset has 85 points with Q 2 < 0.04 GeV 2 . 

We will fit to three types of parameterization. The first is a simple Taylor series expansion, 

2 / 2 \ 2 

q f q 



<W) = l + °ir- + <* (t- +••■> W 

''cut V^cut / 

where we choose to work in units t cut = -im 2 . The second is a continued fraction expansion 
put forward in 



G * {q2) = i ; , = 1 - -n + {aia2 + a ' ] (fS + --- (5) 

1 + Ol -TTTt cut V'cut/ 

We are not aware of a motivation for this ansatz from first principles, but it has been used to 
obtain one of the widely quoted values of the proton charge radius from electron scattering. 
The third is the z expansion described in the Introduction, 

GW) = , + ^ + M + ... = 1 . ^£ + (_| + (£) 2 + ... , (6) 

where z(q 2 ) = z(q 2 ,t cut ,to = 0). As explained below, the coefficients in this expansion are 
bounded; for definiteness here we take \a,k\ < 10. 
We perform fits by minimizing a x 2 function, 

X 2 = ^^(dataj — theory i )i?^ 1 (dataj — theory j) , (7) 

id 

where the error matrix is formed by adding in quadrature the quoted statistical errors, assumed 
uncorrelated, and normalization error, assumed fully correlated within each dataset. In the 
notation of Table 1 of Ref. [19] we use for each experiment, (note that 5 norm refers to the error 
in the cross section) 

\2 



Eij = (5G E )iO~ij + (5 n orm/2) (G E )i(G E ) 3 



4 We obtain similar results by floating the normalization of each experiment and constraining the scale 
factors by an additional contribution to % 2 (as done in |19j ) or by performing the fits at fixed (unit) normal- 
ization and assigning an additional error obtained by adding in quadrature the shift induced by redoing the 
fits with shifted normalization (as done in |20)V 
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polynomial 
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95911 


1 1 99+ 122 




X 2 = 34.49 


32.51 


32.51 


31.10 


28.99 


continued fraction 


882+\° 
X 2 =32.81 


869l| 
32.51 
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z expansion (no bound) 


918±| 


86811 


879l£ 


1022l}° 2 


1193l^ 2 




X 2 =36.14 


32.52 


32.48 


30.35 


28.92 


z expansion ( < 10) 


9181° 


86811 


87911 


880!^ 


8801^ 




X 2 =36.14 


32.52 


32.48 


32.46 


32.45 



Table 1: Proton charge radius extracted from data of Table 1 of [19] (Q 2 < 0.04 GeV 2 ) in 
units of 10~ 18 m, using different functional behaviors of the form factor. Dashes denote fits 
that do not constrain the slope to be positive. 

Errors for the form factor slope are computed by finding the A% 2 = 1 ranged). 

As can be seen from Table [lj the fits with one free parameter differ by many standard 
deviations. Fits with two free parameters agree well, while fits with three or more parameters 
become increasingly unconstrained for the polynomial and continued fraction expansions, as 
well as for the z expansion when no constraints on the expansion coefficients are in place. 
In particular, for k max > 3 in the continued fraction expansion, no meaningful fit can be 
performed (e.g., the slope is not constrained to be positive). 

These results illustrate the problem to be addressed: without detailed knowledge of the 
functional behavior of the form factor, we risk using either too few parameters and biasing the 
fit; or too many parameters and losing predictive power. Note that performing trial fits on 
model data as in [20] is also problematic; some assumption must be made on the functional 
behavior of the form factor in creating the model datasets. To make model independent 
statements requires identifying a bounded class of functions that is guaranteed to contain the 
true form factor, yet is sufficiently restrictive to retain predictive power. The following section 
describes such a class of functions. 

3 Dispersive bounds 

The above fit to the z expansion with a bound on the coefficients illustrates our basic method- 
ology. The present section justifies the \at\ < 10 bound, and demonstrates how further 
constraints can be obtained by disentangling the isoscalar and isovector components of the 
form factor. 

5 We have performed these computations in both MAPLE and MATHEMATICA , and have also checked 
our results using MINOS errors in MINUET. 
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3.1 Form factor definitions 



For completeness we list definitions of the various form factors. The Dirac and Pauli form 
factors, Ff and , respectively, are defined by [211 



(N(p')\J™\N(p))=u(p') 



2m 



^F 2 N (q*) q » 

N 



u(p) 



(9) 



where q 2 = (p' — p) 2 = t and N stands for p or n. The Sachs electric and magnetic form 
factors are related to the Dirac-Pauli basis by [23] 

t 



G N E {t)=F[\t) 



Am 



2 r 2 

N 



N 



At t = they are G E (0) = 1, G n E (0) = 0, G P M (Q) 
write the isoscalar and isovector form factors as 



Fi(t) + F 2 N (t) . 
2.793, G n M (0) = fx r 



(10) 
1.913. We 



G 



(0) 



^ E "r < - J ' 



G 



(1) 



r~<P r 1 



E ■ 



such that at t = they are, G ( E \o) = 1, G%(0) = 1, G^(0) = // p + //„ G^ ; (0) = // p 



»(o), 



Notice that G 



(o) 



(11) 



^ E,Mi ^E,M 



2G V EM for G%\, of [26]. 



3.2 Dispersive bounds 

The analytic structure in the t plane illustrated in Fig. [1] implies the dispersion relation, 

flj^ir ^^") , (12) 

^ tent 

Knowledge of ImG^. over the cut translates into information about the coefficients in the z 
expansion. We begin with a general discussion of these relations. 

Let us consider a general function with the analytic structure as in Fig.[Q G(t) = J2T=o a k z (t) k . 
Equation (j5J) maps points just above (below) the cut in the t plane onto points in the lower 
(upper) half unit circle in the z plane. Parameterizing the unit circle by z(t) = e %e and solving 
(E} for t, we find 

t = t + 2 f cut ~ y = t(9) . (13) 
1 — cos u 

We can now use the orthogonality of z k over the unit circle to find 

a k = -[ d6ReG[t(6) + zO] cos(k6) - - f dQImG[t(6) + zO] sm(k6) . (14) 
Jo 7T J 

Since G is analytic, a k = for k < 0, and therefore 

a = - [ d6Re G[t(6) + zO] = G(t ) , 
7T Jo 

Of" 9 r// It — f 

a k = / dO]mG[t{0) + iQ]sm(ke) = - —— J— — -ImG(t) sm[k6(t)} , fc > 1 . 

7T J0 ^ itcut t-to\ t~ tent 

(15) 
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W\ 

E I 



2/G^(to) 

2 /G f E ) (t ) 



i = t = C t (°- 5GeV 



|G 



^OPE 



||0«G«|| 2 M 1 )(t o )G«(to) 



7.6 
2.5 
14.4 
4.6 



12.1 
3.9 
23.5 
6.7 



tout 


-to 




"tcut 


tcut 


-to 



7=0 



1=1 



1.3 
0.78 



1.8 
1.3 



Table 2: Typical bounds on the coefficient ratios \f^2 k a fc/ a o (upper part of table) and | /ao | 
(lower part) in a vector dominance ansatz. 0ope is defined in Eq. (!23|) . 

The coefficients in the expansion ([3]) can also be used to construct a norm of the form 
factor in the mathematical sense. To keep the discussion general, let us introduce a function 
(f) sharing the domain of analyticity of G, and write 

oo 

<PG = J2 a kZ k . (16) 
Consider the class of norms specified by 

\m\ P = (Eh p Y • ( 17 ) 

In particular, the "uniform norm" is equal to the maximum coefficient size, | |0G| |oo = sup fc |a^| = 
linip^oo ||0G|| p . The case p = 2 is of special interest since the norm is easily related to a dis- 
persion integral, 

The finiteness of | \4>G\ | 2 shows that the coefficients a>k are not only bounded, but must decrease 
in size for sufficiently large k. The relation ||</>G||oo < ||</>G||2 indicates that ||0G|| 2 may 
overestimate the actual size of the relevant coefficients in certain cases. We proceed to consider 
a vector dominance model to illustrate this feature and then turn to a more detailed analysis 
of the spectral functions. 

3.3 Vector dominance ansatz 

In many applications, the || • || 2 norm is used in conjunction with "unitarity bounds" obtained 
by identifying the dispersive integral with a physical production rate. In the present example, 
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dominant contributions to the integral arise from the region below the two-nucleon production 
threshold, and we must turn to different methods of analysis. For example, employing a vector 
dominance ansatz in the appropriate channel, Table [2] displays estimates for the quantity 
| \4>G\ \i/(j>(tn)G(to) = y/J2k a k/ a oi f° r different choices of the functional form of <fi and the 
value of to [j. The effects of the leading resonance in each channel are represented by a Breit 
Wigner profile 



F ( J =°) ^ ^ (19) 

ml-t- iT^m^ m \-t- il p m p 



with oti ~ 1, a 2 ~ —0.12, m w = 783 MeV, r w = 8.5 MeV for the isoscalar channel; and j3i ~ 1, 
fii ~ 3.7, m p = 775 MeV, T p = 149 MeV for the isovector channel. At T = 0, the ansatz is 
normalized to the t = values in Section 13. 11 

We note that in the isoscalar case, the rather large size of the estimated norm is due to the 
narrow width of the u resonance; in fact, in the limit of an infinitely narrow resonance, the 
quantity \\G\\2 diverges, as seen from ()18p . Closer examination indicates that the large norm 
is due not to the coefficients growing in size, but rather to a sequence of coefficients whose slow 
fall-off causes a slow convergence for the sum Y2k a k- A straightforward computation shows 
that the expansion coefficients for an infinitely narrow pole, G(t) = G(0)/ (1 — t/m v ), are for 
k > 1, 



a 




sin 



2k arcsin 




(20) 



In particular, |afc/ao| < 1\J (t cut — 1 ) / (rriy — t cut ). This approximation to the uniform norm 
is also displayed in Table EJ 

Equations ( fl~5l) and ( 1T8]) are model- independent, whereas the approximations based on the 
vector dominance ansatz employed in Table [2] are model dependent. This ansatz aims simply 
to capture the order of magnitude of the coefficients, which is sufficient in practice to constrain 
the form factor fits. The conclusion is that ja^l < 10 is a very conservative estimate for this 
ansatz. 



3.4 Explicit 7T7T continuum 

We can be more explicit in the case of the isovector form factor expansion, where the lead- 
ing singularities are due to hit continuum contributions that are in principle constrained by 
measured irir production and irn — > NN annihilation rates [61 [251 12S] : 

ImGgty) = (t/4 - mlf*F n (tyfl(t) , (21) 

where F n (t) is the pion form factor (normalized according to F n (0) = 1) and is a partial 
amplitude for nn — > NN. Using that these quantities share the same phase we may 
substitute absolute values. Strictly speaking, this relation holds up to the four-pion threshold, 

6 For this purpose we estimate G^io) using a dipolc ansatz for the form factor, G^;(t) ~ 1/(1 — 
t/0.71GeV 2 ) 2 . 
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t < 16m 2 . For the purposes of estimating coefficient bounds, we will take the extension 
of (121 p assuming phase equality through the p peak as a model for the total tttt continuum 
contribution. 

For |-F,r(£)| we take an interpolation using the four t values close to production threshold 
from [22] (0.101 to 0.178 GeV 2 ), and 43 t values from [28] (0.185 to 0.94 GeV 2 ). Values for 
f+(t) are taken from Table 2.4.6.1 of [29]. Evaluating f JT5|) using (12 ip and the experimental 
data up to t = 0.8 GeV 2 ~ 40 m 2 yields for the first few coefficients, at = 1 and to — 0: 
a ~ 2.1 a\ « —1.4, 02 ~ —1.6, 03 ~ —0.9, « 0.2. Using | sm(k8)\ < 1 in the integral gives 
\a k \ < 2.0 for k > 1. 

The leading singularities in the isoscalar channel could in principle be analyzed using data 
for the 37r continuum. Since we do not attempt to raise the isoscalar threshold in our analysis, 
we content ourselves with a simple vector dominance model to estimate the coefficient bounds. 
The first few coefficients for the isoscalar form factor using (1201) for a narrow u resonance are: 
ao = 1, ai ~ —1.2, a 2 ~ —0.96, a 3 « 0.4, 04 ~ 1.3. We will compare the above values to those 
extracted from electron scattering data later. For the moment we note that a bound \a,k\ < 10 
is conservative. 



3.5 Choice of 

Let us return to the choice of 0. We will consider three essentially different choices. First, 
<p(t) = 1 is our default choice. We noted that for = 1 the dominant contributions to ||0G|| 2 
are from narrow resonances. We could negate the large contribution of the leading resonances 
by using for the inverse of a vector meson dominance (VMD) form factor. As a second 
choice, consider 



'VMD 



(t) 



m 



v 



t)/m 



2 

V I 



(22) 



where my is the mass of the leading resonance in the appropriate channel, i.e., p(770) for 
the isovector, w(780) for the isoscalar. Note that using Ge ~ 1/t 2 at large t, the dispersion 
integral remains convergent. There is no loss of model-independence here, since corrections to 
vector dominance are accounted for in the coefficients a^. As discussed in Section [3TB1 a third 
choice of <p is motivated by unitarity and an operator product expansion (OPE): 



'OPE 



(t) 



m N (t cut - t) 4 



6?r (t cut - to)* 



z(t,t mt ,0) 



z(t, teat j to) 
tn-t 



z(t, t cu t, Q 



2 ) 

OPE) 



-Q 



2 

OPE 



(Am 2 N 



0s 

(23) 



where t cut is appropriate to the chosen isospin channel. For definiteness, we choose Qqpe 



1 GeV in the unitarity-inspired 0. In our final fits, we focus on <p = 1 and tc 
demonstrate that the results are essentially unchanged for different choices. 



but 



3.6 Bounds on the region t > 4m 2 



N 

2 



The contribution of the physical region t > Am N to ||</>Ge||2 is 

dt I t cu t — to 



6\\cj)G 



E\\2 



1 



ft J Ami, t ~ ^0 V t — t cn t 



\<1>G E \ 2 



(24) 



The cross section for e + e — » iViV is [30 



a(f) 



47ra 2 
~3T 



1 - 



4m|r 



Gm(*)| j 



and thus for the proton electric form factor we have 



t 



s\\^ E \\l 



1 



7T 



dt 



inn 



cut 



to 



t — to V t — t. 



a(t) 



[a (t)v(t) \G M /G! 



2m 2 N /t 



(25) 



(26) 



where ao = Ana 2 /?>t and v(t) = yl - Am? N /t is the nucleon velocity in the center-of-mass 
frame. Using the data from [31] (see also [32U33]). we can perform the integral from t = 
4.0 GeV 2 to 9.4 GeV 2 assuming \G P M /G P E \ < l@At t = and <p = 1, we find the result 
^II^jeII! ~ (0.03) 2 , to be added to the contribution from t < Am 2 N . This result is obtained 
by using for a(t) the measured central value plus la error. The remaining integral above 
t = 9.4 GeV 2 can be conservatively estimated by assuming a constant form factor beyond this 
point, yielding an additional 5||G^||| ~ (0.008) 2 . The neutron form factor can be treated 
similarly using the data from [M] for t = 3.61 to 5.95 GeV 2 . This leads to 5\\G n E \\ 2 2 w (0.05) 2 . 
The remainder at high t assuming a constant form factor yields an additional ^HG^HI ~ 
(0.05) 2 . Similarly, using llmG^ sin k0\ < \Ge\ the contribution of the timelike region to (I15p 
is small: \8a k \ < 0.011 + 0.004 for the proton, and \8a k \ < 0.013 + 0.025 for the neutron. We 
conclude that when estimating the bounds on coefficients, the physical timelike region can be 
safely neglected. 

Let us mention that we can bound the contribution of the physical timelike region by 
a perturbative quark-level computation. Decompose the electromagnetic current correlation 
function as 



d*xe^(0\T{J^(x) t -C(0)}|0> = (<f <f - tfgrmq 



and define 



x(Q 



l d 2 



OPE) 



2 8{q 



2\2 



(? 2 n(g 2 )) 



1 f°° , tlmU(t) 
dt 



q 2 =~Q 



2 

OPE 



7T 



Jo 



(t + Q 



(27) 



(28) 



OPE J 



The two-nucleon contribution to the correlator satisfies 



ImIT(t) > 



in 



N 



Gnt 



1-^\<PG E \ 2 



(29) 



7 For \Gm/Ge \ > 1- the quantity in square brackets in (f26| is bounded by the quantity denoted by |G| 2 in 
[31] . This inequality is satisfied experimentally in the t range of interest. 
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and hence with 4>Ge = J2 k a>kZ k and the choice of <p in (1231) . 



x(Qope) > 



TT 



1 




'JV 



|0G E | 2 > 8\\<j)G E 



(30) 



If we choose Qq PE l &r ge enough, the function x(Qope) i s perturbatively calculable as an 
operator product expansion: x ~ X]/ e //8 7r2( 3oPE a ^ leading order, where e/ denotes the 
electric charge of a given quark flavor. Choosing for illustration Qq PE = 1 GeV 2 , nj = 

3 light quark flavors, and i cut = 4m 2 , we find the bounds 8(J2k a t) ~ (1-0) 2 f° r = 
and 5(Efc«fc/«o) ~ ( x - 4 ) 2 for *o = C*( - 5 GeV 2 ). We note that these "unitarity bounds" 
overestimate the contribution from the physical region t > 4:m 2 N , due both to subthreshold 
resonance production, and to other channels, e.g., NN plus pions, above threshold. For this 
reason, we do not dwell on a more precise analysis of this bound, or on a separation into 
definite isospin channels. 

4 Proton charge radius extraction 

We consider several possibilities to reduce the error bars for the proton charge radius ex- 
tracted in Section [2J We first consider the inclusion of higher-Q 2 data. We then optimize the 
charge radius extraction by separating isoscalar and isovector components, recognizing that 
the isoscalar threshold is at 9m 2 . At the same time, we illustrate the (small) effect of differ- 
ent expansion schemes. Finally, we consider the possibility to effectively raise the isovector 
threshold by constraining the spectral function between 4m 2 and 16m 2 . 

4.1 Including higher Q 2 data 

We have argued that, taking the data tabulated in p2] at face value, the final entry in Tabled] 
is a model-independent determination of the proton charge radius: r v E = 0.878^o 062 ^ n 
the absence of further model-independent constraints on the form factors, obtaining a proton 
charge radius with smaller error requires further experimental input. Here we investigate the 
impact of higher-Q 2 proton scattering data. 

Figure [3] shows the central value and la (Ax 2 = 1) error band obtained by fitting the 
electron-proton scattering data compiled by Arrington et al. [7]. We take = 1 and to = 0, and 
include as many coefficients aj, as necessary for the fits to stabilize. As the figure illustrates, 
for Q 2 > few x 0.1 GeV 2 the impact of additional data is minimal. While an ever greater 
number of coefficients a k at higher k must be included to obtain convergence, the total error 
on the slope at Q 2 = is not reduced. For later use, we note that the coefficients 0^=1,2,3 
extracted from the fit at Ql mx = 1 GeV 2 are -1.01(6), -1.41$, 2^ 2 . 

4.2 Raising the isoscalar threshold: inclusion of neutron data 

We can separate the isoscalar from the isovector form factor, making use of the fact that the 
isoscalar cut is further away from t = than the isovector cut, translating to a smaller value 
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Figure 3: Variation of the fitted proton charge radius as a function of maximum Q 2 . Fits of 
the proton data were performed with k ma , x = 10, <p = 1, to = 0, |afc| < 10. Data from [7J. 



of Mmax as discussed in the Introduction. A combined fit of proton and neutron data can then 
be performed. For the proton form factor we again use the data from [7J. For the neutron 
electric form factor, we use 20 data points from [35l [36l EZl [38l EHl SQl SU S21 SHI HH HSl S6] . 
We take as additional input the neutron charge radius from neutron-electron scattering length 
measurements [5]: 

(r 2 )™ = -0.1161(22) fm 2 . (31) 

Table [3] shows the effect of different expansion schemes (choices of and to) an d coefficient 
bounds on the form factor slope determination. For later use, the coefficients afc=i,2,3 extracted 
from the fit for Q 2 max = 1 GeV 2 , = 1, t = and k max = 8 are -1.99^;^, 0.3+};jj, -2ljj for 
the isoscalar channel; and — 1.20ioo5> — 0.6Zt^|, — 2ly for the isovector channel. The sign and 
approximate magnitude of the first coefficients agree with the tttt continuum model, and the 
narrow- width u resonance model mentioned in Section 13.41 

4.3 Raising the isovector threshold: inclusion of jnr data 

We can effectively raise the isovector threshold by including the tttt continuum explicitly, as 
constrained by tttt production and hit — > NN data: 

G { i\t) = G cut (t) + ^2a k z k (t,t cnt = 16m 2 , t ) , (32) 

k 

where G cut (t) is generated by (]2~T|) for 4m 2 < t < 16m 2 . For 1^(^)1 we take the four t 
values close to production threshold from [27] (0.101 to 0.178 GeV 2 ), and twelve t values 
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k - 9 


3 


4 


5 


6 


= 1, to = 0, 


\a k \ < 10 


888+J 


8651U 


888I22 




878l 2 ° 






X 2 = 33.67 


23.65 


21.80 


21.13 


20.47 


= 1, t = 0, 


\a k \ < 5 


8881J 


865+jl 


88111° 


8851^ 


55Z_2o 






X 2 = 33.67 


23.65 


21.95 


21.46 


21.06 


= 0VMD; ^0 


= 0, \a k \ < 10 


865±| 


874l^ 2 


884l| 


879+1 


8771 22 






X 2 = 23.26 


22.50 


22.15 


21.59 


21.09 


= 1, t = 




888lj? 


865±\\ 


8801^ 


88211* 


882±1| 






X 2 = 33.67 


23.65 


22.07 


21.45 


21.18 


= 0OPE, t = 


= 


904+J 


86111° 


888!^ 


qqq+20 
88d_ 20 


881! 20 , 






X 2 = 61.34 


24.38 


21.62 


20.86 


20.51 


= 0OPE, *0 = 


= t° pt (0.5GeV 2 ) 


912±| 


869l| 


8871H 


88ll 2 ° 


8801 20 , 






X 2 = 93.69 


22.54 


21.05 


20.32 


20.32 



Table 3: RMS charge radius extracted using electron-proton and electron- neutron scattering 
data, and different schemes presented in the text. The neutron form factor slope is constrained 
using (|3T|) . A cut Q 2 ^^ = 0.5 GeV 2 is enforced. In the lower part of the table, the bounds on 
J2 k al from Table [2] are multiplied by 4. 0vmd and 0ope are defined in Eqs. (l2~2i) . (l2~3l . 

from [28] (0.185 to 0.314 GeV 2 ). The product of the remaining kinematic factor and f} 
from [29] is interpolated to the appropriate t value, and the integral computed as a discrete 
sum. Using coarser bin size (e.g. 8 instead of 16 bins) has no significant effect, indicating 
that discretization error is small. Estimating the remaining coefficients by modeling the tttt 
continuum contribution for 16m 2 < t < 40m 2 using (fT5l) and ([21]) at = 1 and to = gives 
coefficients a\ ~ —4.5, a 2 ~ 2.2, 0,3 ~ 2.1. Setting |sin(/c#)| in ([15]) yields \a k \ < 5.0 for the 
remaining contribution of the nir continuum in this model. 

We fit using the same proton and neutron data as in Section 14.21 The resulting fit coeffi- 
cients a fe= i |2>3 for Q 2 max = 1 GeV 2 , = 1, t = and k max = 8 are -1.93(6), -0.5l^, 2 ± 7 for 
the isocalar form factor; and — 3.401q'5q, 3.7l}"3, 3-io for the isovector form factor. The sign 
and approximate magnitude of the first coefficients agree with the remaining nit continuum 
model discussed above in the isovector case; and with the oj pole model discussed at the end 
of Section 13.41 for the isoscalar case. The sizable contribution of the isovector a k =i in this 
scheme can be traced to the residual effects of the tttt continuum, including the p peak, near 
the higher threshold. With no loss of model-independence, we can replace G cn t(t) above with 
a new G C ut(^) generated by (12 ip for 4m 2 < t < 40m 2 , i.e., with the hit continuum modeled 
to larger t. The value t cut = 16m 2 remains the same. We emphasize that this does not in- 
troduce a model dependence, as any discrepancy between G cut (t) and the true mi continuum 
is accounted for by parameters in the z expansion. The resulting central value and errors on 
the charge radius are changed minimally by this modification. The isoscalar coefficients are 
also not significantly changed, while the isovector coefficients become 1.07(10), 1.6ll;5, l-s- 
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Figure 4: Variation of the fitted proton charge radius as a function of maximum Q 2 . Fits 
were performed including proton data, neutron data and the mr continuum contribution to 
the isovector spectral function, as detailed in the text. Fits were performed with /c max = 8, 
0=1, t = 0, \a k \ <10. 



Figure H] shows the resulting extraction of the proton charge radius using for G cut (t) the full 
model of the irn continuum, and our default <p — 1, to — 0. As in Fig. [3J the inclusion of data 
beyond Q 2 ~ few x 0.1 GeV 2 has minimal impact on the fits. 



5 Discussion 

We have discussed determinations of the proton charge radius from the slope of the proton 
form factor G p E (t), in four cases: (1) low-Q 2 electron-proton scattering data; (2) proton data 
including high Q 2 \ (3) proton plus neutron data; and (4) proton, neutron, and im data. We 
have investigated various expansion schemes, corresponding to choices of the parameter t and 
the function 0, and shown that the impact on r v E is minimal; in the following discussion we 
take = 1 and t = 0. 

Including just the low Q 2 proton data [19], we find the result as in Table [1] [case (1)] 
r E = 0.877^q'q49 ± 0.011 fm, where the first error is obtained using the more stringent bound 
|afe| < 5, and the additional error is conservatively estimated by finding the maximum variation 
of the Ax 2 = 1 interval when the fits are redone assuming \a,k\ < 10. Using a larger Q 2 range 
of proton data [7] decreases the uncertainty. Taking for definiteness Q max = 0.5 GeV 2 and 
fc max = 8, we obtain via the same procedure, as in Fig. [3] [case (2)] r p E = 0.870±0.023±0.012 fm. 
Including the neutron data, as in Table El we find [case (3)] r v E = 0. 880^020 ± 0.007 fm, 
where the same bounds, ja^l < 5, \a k \ < 10 are enforced on both isoscalar and isovector 
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coefficients and again k mSLX = 8^ The uncertainty induced by the neutron charge radius (13~T|) 
is negligible in comparison, < 0.0005 fm. Finally, including G cut (t) as in (I3"2"j) . we find [case 
(4)] r p E = 0.871 ± 0.009 ± 0.002 ± 0.002 fm. For definiteness, we here include in G cut (t) the 
extension of the nir continuum model up to t = 40 m 2 . The first and the second error are as 
above, and the final error is obtained by assigning a 30% normalization error to the continuum 
contribution, as discussed below. 

Let us compare our results to several previous determinations of r p E . Many of these suffer 
from model assumptions on the functional behavior of the form factor. The small uncertainties 
obtained by Simon et al. 07] (r p E = 0.862 ±0.012) and by Rosenfelder [19] (r p E = 0.880 ±0.015) 
require inputs from higher Q 2 data, which however we do not believe were robustly estimated. 
We find that the coefficient of t 2 in the expansion of G p E (t) [Eq.(j4])] is constrained by the 
Arrington et al. data compilation [7] to be aJ aylor /^ ut = 0.014^;^ ± 0.005 fm 4 (using Q max = 
lGeV 2 , & max = 10). A much smaller uncertainty, aJ aylor /t 2 cut = 0.011(4) fm 4 or 0.014(4) fm 4 , 
was adopted in [19]. Even neglecting the additional uncertainty due to cubic and higher 
order terms, this would lead to a result 0.878 ± 0.0081o;q 4 q obtained using (0J and data as 
in Table [TJ The errors are from the data and from the first uncertainty on the quadratic 
coefficient, respectively. 

The analyses of Sick [20] (r p E = 0.895 ± 0.010 ± 0.013) and Blunden and Sick [H] (r p E = 
0.897 ± 0.018) employ the continued fraction expansion ([5]). This functional form is unstable 
to the inclusion of additional parameters (see Table [TJ), and error estimation relies on the 
investigation of model datasets. In this paper we have not fit directly to cross section data, 
and we have not applied our analysis to this dataset. For a variation of this analysis see [19]. 

The dispersion analysis of Belushkin et al. [26] (r p E = 0.8441°;°°* fm, r p E = 0.830ig;gg|) 
does not attempt to estimate uncertainties due to the constrained shape of the assumed form 
factors. Our analysis makes clear which inputs have the most effect on the charge radius 
extractions. In particular, data at large for either timelike or spacelike t, has minimal 
impact on fits to obtain Q 2 « quantities. Inclusion of high-Q 2 data does introduce sensitivity 
to additional parameters, whose omission would introduce model dependence. Our analysis 
provides a systematic procedure to analyze a wide range of datasets in a model-independent 
way. We emphasize that our goal is not simply reduction in the quoted error, but also the 
robust estimation of uncertainties. 

Regarding the bounds on coefficients, in all approximations that we have considered the 
bound |afe| < 10 appears very conservative. The sign and magnitudes of the first coefficients 
are consistent with expectations based on simple models, and it is rigorously true that the 
coefficients at must eventually decrease in magnitude for large k. At a practical level, the 
experimental determinations of these coefficients in each of the cases (l)-(4) above are consis- 
tent with magnitudes not larger than |eifc| ~ 2. Our implementation of the bounds on could 
be formalized in terms of standard methods of constrained curve fitting [50J . As discussed in 
[T5] , our assumption of a flat "prior" should be conservative. 

Our analysis cannot discern inaccuracies in the datasets. For example, we have assumed 
that radiative corrections are properly accounted for in the compilations [191 U]> an d that 

8 The slight difference between this value and that inferred from the final column for the first two rows of 
Table [3] is due to the slight difference between fc max = 6 and fc max = 8. 
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data correlations are sufficiently described by our treatment Within these assumptions, the 
values for cases (l)-(3) represent model-independent determinations of the form-factor slope. 
Case (4) is more subtle. While ( 12 ip is a model-independent relation for the stated range 
4m 2 < t < 16m 2 , the determination of f]_(t) in this range involves a dispersion relation with 
contributions from values of t where the function is not rigorously constrained by continuation 
of nN scattering data. Errors are not given in the tabulation [29], and we are not aware of 
a critical assessment of the uncertainties associated with this analysis. It may be interesting 
to revisit this question. Ref. [21] suggests a 15% error in the normalization of /+(£) at the p 
peak; we take twice this value, 30%, as a representative uncertainty, which encompasses also 
the errors in The resulting error for case (4) is thus not as rigorous, although the 

resulting /+(£) would need to be very different to become a dominant source of error. 

Within the stated uncertainties we find consistent results in each of our determinations, 
using both low and high Q 2 proton data, neutron data, and pion continuum data. These 
methods can be applied to other datasets, and to fits using partial cross sections versus ex- 
tracted form factors. For example, in a recent set of results [53] the variation of r p E under 
different model shapes for the form factor is larger than the other stated statistical and sys- 
tematic errors. The same methods can be applied to other nucleon form factors and derived 
observables, including the axial- vector form factor probed in neutrino scattering |54j . 
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